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An  interpretation  scheme  has  been  developed  for  the  determination  of 
laterally  homogeneous  earth  structure  using  body-wave  travel  time  and  ampli¬ 
tude  data.  A  family  of  p { x )  (ray  parameter  v^  distance)  functions  consistent 
with  the  data  are  obtained  as  solutions  to  an  inverse  problem.  Each  p(x) 
solution  has  a  corresponding  velocity-depth  function  which  maximizes  or 
minimizes  the  depth  for  one  velocity.  The  family  of  allowable  p(x)  curves 
therefore  defines  a  velocity-depth  envelope  which  bounds  the  range  of  earth 
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SUMMARY 

An  interpretation  technique  has  been  developed  for  the  determina¬ 
tion  of  laterally  homogeneous  earth  structure  using  body-wave  travel  time 
and  amplitude  data.  A  family  of  p(x)  (ray  parameter  vs^  distance) 
functions  consistent  with  the  data  are  obtained  as  solutions  to  an  inverse 
problem.  Each  p(x)  solution  has  a  corresponding  velocity-depth  function 
which  maximizes  or  minimizes  the  depth  for  one  velocity.  The  family  of 
allowable  p(x)  curves  therefore  defines  a  velocity-depth  envelope 
which  bounds  the  range  of  earth  structures  consistent  with  these  data. 

1  INTRODUCTION 

Earth  structure  cannot  be  determined  uniquely  from  a  finite  set 
of  body-wave  observations  containing  errors.  In  general,  there  is  a 
broad  range  of  velocity-depth  models  which  are  consistent  with  a  given 
data  set.  In  recent  years,  several  authors,  most  notably  Bessonova,  et  al . 
(1974,  1976),  have  proposed  methods  for  estimating  bounds  for  the  allowable 
earth  structures  by  inverting  travel  time  data.  However,  since  these 
techniques  neglect  amplitude  information,  the  resulting  bounds  are  somewhat 
less  restrictive  than  may  be  indicated  by  the  observed  data.  Mellman  (1978) 
has  taken  a  different  approach,  whereby  velocity-depth  models  are  deter¬ 
mined  by  matching  the  waveforms  of  observed  and  synthetic  seismograms. 

The  non-uniqueness  of  the  resulting  earth  models  is  estimated  from  the 
discrepancy  between  solutions  obtained  by  using  different  starting  models 
in  the  inversion  scheme. 
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In  the  following  conmunication,  we  propose  a  technique  for  inter¬ 
preting  seismic  refraction  data  which  incorporates  both  amplitude  and 
travel  time  information.  The  data  are  inverted  to  determine  the  range  of 
seismic  velocity-depth  models  which  are  consistent  with  these  observations. 
The  earth  models  in  this  investigation  are  restricted  to  be  laterally 
uniform,  since  we  allow  the  velocity  to  vary  only  as  a  function  of  depth. 
Furthermore,  this  investigation  is  restricted  to  velocity-depth  models 
which  contain  no  low-velocity  zones. 

The  central  idea  in  this  interpretation  scheme  is  the  choice  of 
the  function  p(x)  (wave  slowness  as  a  function  of  source-receiver 
distance)  for  performing  the  inversion.  There  are  several  advantages 


for  deforming  the  p(x)  curve  to  fit  the  data: 

(a)  The  travel  time  at  any  distance  is  a  linear  integral 


a  cl..- _ ; 

I<  j  ii  . 


of  the  p(x)  curve.  , 

(b)  The  velocity-depth  (v(z))  model  can  be  obtained  by  i,  . 

a  simple  linear  Weichert-Herglotz  integral  of  the  -r:'‘  1 

A  . .  : 

p(x)  curve.  i , 

T-.„.  j‘‘ 

(c)  Amplitude  information  can  be  recovered  from  the  p(x)  » 
curve  through  a  simple  Disc  Ray  Theory  integration  1  j 
(Wiggins,  1976). 

Hence,  aVl^  of  the  desired  transformations  can  be  made  directly  from  the 
p(x)  curve.  These  transformations  are  well-defined  and  stable  functions 
of  the  p(x)  curve. 
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Another  Important  reason  for  choosing  the  p(x)  function  for  per¬ 
forming  the  inversion  is  based  on  the  need  for  a  simple  representation 
of  the  non-uniqueness.  McMechan  &  Wiggins  (1972)  and  Wiggins,  et  al. 
(1973)  have  shown  how  the  uncertainty  of  the  earth  structure  can  be 
represented  by  an  envelope  in  the  p-x  plane,  which  bounds  the  set  of 
allowable  p(x)  curves.  They  have  also  shown  how  the  p(x)  function 
can  be  manipulated  in  order  to  determine  extremal  depths  corresponding 
to  any  velocity.  Their  methodology  will  be  used  in  this  interpretation 
scheme  and  will,  in  fact,  be  incorporated  directly  into  the  inverse 
problem  so  that  extremal  depths  can  be  determined  for  the  set  of  per- 
missable  velocity  models. 

The  interpretation  procedure  involves  a  combination  of  forward  and 
inverse  modeling.  The  inverse  problem  consists  of  determining  the  set  of 
p(x)  curves  for  which  the  velocity-depth  functions  are  extremal.  The 
forward  problem  consists  of  determining  the  synthetic  seismograms  corres¬ 
ponding  to  these  extremal  p(x)  solutions.  In  general,  all  of  the  data 
cannot  be  fit  in  one  iteration,  so  the  problem  proceeds  iteratively.  The 
seismograms  corresponding  to  the  initial  p(x)  solutions  have  amplitudes 
and  travel  times  which  are  inconsistent  with  the  observed  data.  This  can 
be  remedied  in  subsequent  iterations  by  imposing  limits  which  restrict 
the  p(x)  curves  from  certain  regions  of  the  p-x  plane.  These  bounds 
form  a  p(x)  envelope  which  is  refined  on  each  iteration  until  all  of 
the  solutions  match  the  observed  data  to  within  some  uncertainty.  The 
scatter  of  the  calculated  seismograms  is  governed  by  the  width  of  the 
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p(x)  envelope,  which,  in  turn,  determines  the  range  of  allowable  earth 
structures.  , 

The  interpretation  scheme  will  be  illustrated  in  detail  in  a 
later  section.  First,  we  turn  our  attention  to  the  most  complex  aspect 
of  the  procedure  —  the  inverse  problem. 

2  THE  INVERSE  PROBLEM 

The  inverse  part  of  the  interpretation  scheme  consists  of  deter¬ 
mining  an  optimal  p(x)  curve  which  satisfies  certain  travel  time 
constraints  and  at  the  same  time  maximizes  or  minimizes  the  depth  for 
some  velocity.  Both  of  these  conditions  are  quite  straightforward  in 
that  they  involve  transformations  from  the  T-x  and  v-z  to  the  p-x 
domains.  The  transformations  are  linear  and  are  written  mathematically 
as  equalities.  However,  there  are  also  a  number  of  bounding  or  inequality 
constraints  which  we  wish  to  impose  on  the  p(x)  solutions.  Two  of  these 
require  the  velocity-depth  function  corresponding  to  any  p(x)  solution 
to  be  physically  meaningful;  that  is,  the  velocity  must  be  single  valued 
in  depth  and  the  depths  for  all  velocities  must  be  positive.  A  third 
inequality  constraint  stems  from  the  forward  problem,  where  bounds  are 
imposed  on  the  solutions  limiting  them  from  certain  regions  of  the  p-x 
plane.  Least-squares  error  techniques  are  not  well  suited  to  problems 
containing  many  inequality  constraints.  We  have  turned,  therefore,  to 
the  realm  of  L^norm  linear  estimation  in  order  to  solve  our  inverse 
problem.  However,  before  we  discuss  our  method  for  solving  it,  let's 
examine  the  problem  more  closely. 
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We  cast 
with  the  p(x) 


the  inverse  problem  in  the  form  of  a  system  of  equations 
curve  as  a  vector  of  unknowns  x,  i.e.. 


A  x  -  d  -v  e  (1 ) 

where  each  row  has  the  form 
M 

2\j  xj  •  dkiek  ■  k*'-N-  (2) 

j=i 

Here  x  is  the  discretized  p(x)  curve,  d  is  the  "data"  vector,  and 
A  is  the  transfer  matrix  which  relates  the  unknowns  x  to  the  "data" 
d  .  The  optimum  solution  for  x  is  determined  by  minimizing  the  error 
vector  e  in  some  way. 

The  unknowns  in  this  problem  are  chosen  to  be  the  discrete  source- 

receiver  distances  x.,  j  =  1,  M  corresponding  to  finite  samplings  p., 

J  J 

j  =  1,  M  of  the  ray  parameter  over  some  interval  of  interest  [Pm^n>  P^^- 
The  distance  x  is  chosen  to  be  the  dependent  variable  because  it  is 
generally  interpreted  as  being  single-valued  in  p.  (The  converse  is 
not  true  if  there  are  any  triplications  in  the  p(x)  curve.)  We  use 
linear  interpolation  between  the  discrete  points  x^  and  Xj+1  in 
order  to  approximate  a  continuous  p(x)  curve.  We  have  found  that 
this  piecewise  linear  representation  has  significant  benefits  over  a 
constant  interpolation  scheme,  since  discretization  effects  are  reduced. 
(For  simplicity,  the  equations  in  this  communication  are  presented  as 
having  constant  interpolation,  i.e.,  as  simple  summations.) 


6 


The  p  interval  depends  on  the  problem  at  hand  and  must  be  deter¬ 
mined  first.  In  general,  p[nax  is  equal  to  (velocity)”^  at  the  surface  and 

is  assumed  to  be  known.  The  lower  limit  p  .  is  estimated  from  the 

min 

travel  time  data  at  large  source-receiver  distances. 


3  SOLVING  THE  INVERSE  PROBLEM 

The  bounding  constraints  which  we  impose  on  the  solution  mani¬ 
fest  as  inequality  equations.  For  example,  the  equations  restricting  the 
p(x)  curve  to  lie  within  some  limits  x*-0  and  xHI  are  simply 


xj°  <  xj  <  xj 1  ,  j  =  1,  M  .  (3) 

One  way  to  ensure  that  such  bounding  constraints  are  always  satisfied  is 
to  insert  the  relevant  equations  directly  into  the  inverse  problem. 
Claerbout  &  Muir  (1973)  have  devised  an  algorithm  for  solving  an  over¬ 
determined  system  of  equations  that  can  be  used  to  include  inequality 
constraints  into  the  problem.  Their  algorithm  optimizes  the  solution 
according  to  the  absolute  value  error  criteria,  or  L^  norm: 


minimum 


(4) 


in  contrast  with  the  L2  norm  which  minimizes  the  sum  of  the  squares  of 
the  errors.  The  term  in  equation  (4)  is  a  special  weighting  function 
which  may  be  asymmetric  about  e^  =  0.  Each  row  in  the  inverse  problem  (1 ) 
is  assigned  two  weighting  coefficients,  one  for  positive  error  and  the 
other  for  negative  error.  If  w*  =  w”  ,  this  particular  equation  is  a 
simple  equality.  However,  if  one  of  the  weights  is  made  substantially 
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larger  than  the  other,  the  positive  and  negative  errors  are  penalized 
differently  and  the  corresponding  equation  becomes  an  inequality.  This 
asymmetric  error  capability  proves  very  useful  in  many  applications  such 
as  this  one,  where  a  large  number  of  inequality  constraints  are  imposed 
simultaneously  on  the  solution. 

In  this,  as  in  any  inverse  problem,  difficulties  can  arise  if  the 
number  of  unknowns  exceeds  the  number  of  observations,  or  if  the  distri¬ 
bution  of  data  is  biased  in  such  a  way  that  certain  parts  of  the  solution 
are  ill-determined.  In  these  situations,  some  kind  of  smoothing  condition 
must  be  applied  in  order  to  ensure  a  unique  solution.  The  smoothness 
condition  used  in  this  application  is  a  linear  depth  vs  velocity  relation. 
(Since  a  linear  interpolation  is  used  for  the  x(p)  curve  between  successive 
points  x(pi)  and  x(pi+1),  the  corresponding  depth  vs  velocity  curve  cannot 
be  linear.  The  smoothness  condition  applies  only  to  the  sample  points  pi , 
i.e.,  z(l/p._i)  -  2 z ( 1  / p  1- )  +  z(l/p.+1)  =  0.)  The  smoothness  condition  is 

applied  by  simply  adding  the  relevant  equations  to  the  system  of  equations  (1) 

■j' 

and  adjusting  the  weighting  function  to~  in  such  a  way  that  the  data 
constraints  always  dominate.  The  necessary  adjustment  is  that  the  ^  for 
the  linear  depth  vs^  velocity  relations  must  be  less  than  the  weighting 
function  for  the  data  divided  by  the  total  number  of  linear  z(l/p) 
constraints. 

In  its  final  form,  the  inverse  problem  contains  equations  corres¬ 
ponding  to  five  different  types  of  constraint:  travel  time  information,  an 
extremal  depth  condition,  a  single  valued  velocity  restriction,  p-x  bounds, 
and  a  smoothness  criterion.  Clearly,  all  of  these  conditions  together  give 
us  an  over-determined  system  of  equations.  The  algorithm  solves  this 
type  of  problem  by  selecting  M  of  the  N  equations  to  be  satisfied 


exactly  (M  is  the  number  of  unknowns).  We  point  out  that  this  feature  of 
the  L.|  solution  does  not  imply  that  the  remaining  N-M  equations  are 
ignored;  the  particular  set  of  M  equations  to  be  satisfied  are  selected 
so  as  to  minimize  the  cumulative  errors  from  the  problem  as  a  whole. 


4  TRAVEL  TIME  CONSTRAINTS 

The  ray  parameter  p  at  some  range  x  is  related  to  the  travel 
time  curve  through  a  derivative; 


P  =  dlM 
p  dx 

The  travel  time  is  therefore  simply  an  area  under  the  p(x)  curve 
x 

T(x)  =  f  p(y)  dy  . 

"'o 


(5) 


(6) 


Since  p  is  not  always  single  valued  in  x,  it  is  more  convenient 
to  conduct  the  integration  in  p  rather  than  x.  Integration  of 
equation  (6)  by  parts  yields 


_max 

T(x)  =  p  •  x(p)  +  j  x(q)  dq 

P(x) 

=  p  •  x(p)  +  t(p)  (7) 

The  delay  time  x(p)  represents  the  time  intercept  at  x  =  0  for  ray 
parameter  p  and  travel  time  T(x).  For  a  discrete  p(x)  curve,  -r(p) 
is  simply  a  summation  of  the  source-receiver  distances  over  some  interval 
of  ray  parameters  [P.,  PTOX]  , 


T  . 
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(8) 


i 

pmax 

=  E  xo  apj 
j=i 

Equation  (8)  is  in  the  same  form  as  one  row  of  the  inverse  problem, 
equation  (2).  Thus,  t(p)  data  can  be  entered  directly  into  the  inverse 
problem  and  fitted  in  one  iteration.  However,  this  is  not  the  case  for 
T(x)  data.  Travel  times  must  be  converted  to  the  appropriate  integral 
form  using  equation  (7)  before  they  can  be  entered  into  the  inverse  problem. 
Note  that  in  equation  (7),  each  T(x)  pair  requires  a  corresponding  p 
value  in  order  to  set  the  lower  limit  of  integration.  Initial  estimates 
for  p  at  each  observation  distance  x  can  be  obtained  by  differentiating 
the  T(x)  curve.  One  way  to  improve  these  p  values  is  by  solving  the 
inverse  problem  iteratively  without  the  extremal  depth  condition  until  the 
p(x)  solution  matches  the  p  limit  used  to  evaluate  equation  (7). 

In  the  inverse  problem,  T(x)  data  provide  simple  integral  constraints 
on  the  p(x)  curve,  just  as  t(p)  data  do.  Integral  constraints  place 
no  restrictions  on  the  individual  distances  x.,  which  means  that  the  travel 

J 

times  for  the  initial  solution  may  not  match  the  observed  T(x)  values.  In 
other  words,  the  velocity-depth  model  is  not  yet  consistent  with  the  travel 
time  data.  One  way  to  remedy  this  situation  is  to  place  bounds  on  the 
p(x)  solution.  These  bounds  are  determined  by  comparing  the  travel  times 
of  the  inverse  solution  with  those  of  the  observed  seismograms.  In  general, 
this  p(x)  envelope  must  be  refined  by  performing  the  inverse  and  forward 
problems  several  times,  until  a  satisfactory  fit  is  obtained.  The  discrepancy 
between  the  observed  and  calculated  seismograms  depends  on  the  width  of 
the  p(x)  envelope  used  to  constrain  the  solution.  Thus,  two  types  of 
constraints  are  needed  to  match  travel  time  data.  The  procedure  for 
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matching  travel  time  data  will  be  illustrated  further  in  the  following  sections 
when  we  present  some  examples. 

If  there  are  any  triplications  in  the  observed  T(x)  data*  it  is 
advantageous  to  use  the  differential  travel  times  between  the  early  and 
later  arrivals  to  constrain  the  p(x)  solutions.  The  differential  travel 
times  provide  integral  constraints  over  smaller  sections  of  the  p(x) 
curve  and  therefore  provide  stronger  constraints  than  do  the  absolute  arrival 
times.  By  a  similar  argument,  we  expect  the  extremal  velocity-depth  limits 
to  reflect  the  number  of  integral  constraints  imposed  —  in  other  words,  the 
number  of  T(x)  data  used. 

5  EXTREMAL  DEPTH  CONSTRAINT 

McMechan  &  Wiggins  (1972)  and  Wiggins,  et  al.  (1973)  have  shown 
how  the  p(x)  function  can  be  manipulated  in  order  to  obtain  maximum  and 
minimum  depths  corresponding  to  any  given  velocity.  In  this  section, 
we  will  illustrate  how  such  an  extremal  depth  condition  can  be  incorporated 
directly  into  the  inverse  problem.  This  will  automatical ly  give  an  extremal 
earth  model  which  is  consistent  with  the  data  each  time  the  inversion  is 
performed. 

The  depth  z(p)  corresponding  to  the  turning  point  of  ray  p  is 
obtained  by  the  Weichert-Herglotz  integral  of  the  p(x)  curve: 


where  the  corresponding  velocity  is  v  =  1/p.  For  a  discrete  p(x) 
curve,  equation  (9)  becomes 
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which  is  simply  a  weighted  summation  of  the  discrete  distances  x.  over 

J 

the  p  interval  [p^,  p  ].  These  weighting  coefficients  are  positive 
and  decrease  monotonically  in  p.  In  order  to  determine  the  p(x)  curve 
which  satisfies  the  data  constraints  and  gives  an  extremal  depth  for  some 
ray  parameter  p*,  say,  one  row  containing  the  Weichert-Herglotz  coefficients 
for  z(p*)  is  added  to  the  inverse  problem,  equation  (1).  Let's  say  we 
wish  to  minimize  the  depth  for  p*.  The  corresponding  depth  value  in  the 
vector  d  is  then  set  to  zero  and  the  asymmetric  weighting  function  w1 
for  this  equation  is  set  to  some  value  less  than  w*  for  the  data  con¬ 
straints  and  larger  than  the  sum  of  the  w*  for  the  linear  v(z)  relations. 
If  the  weighting  function  cT  is  chosen  correctly,  the  L1  algorithm  will 
find  a  solution  which  accomodates  for  this  extremal  depth  constraint,  but 
does  not  actually  solve  it.  The  data  equations  will  take  precedence. 

Consider  the  following  example. 

6  MONOTONIC  p(x)  EXAMPLE 

Fig.  1(a)  shows  three  seismograms  which  were  computed  from  the 
p(x)  curve  in  Fig.  1(b)  using  Disc  Ray  Theory  integration.  The  solid  curve 
in  Fig.  1(c)  shows  the  corresponding  earth  model.  These  three  seismograms 
provide  three  travel  time  data  which  were  used  as  input  data  for  the  inverse 
problem.  The  travel  times,  along  with  their  corresponding  p  values,  are 
listed  in  Table  1 . 
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To  fix  ideas,  we  first  illustrate  the  set-up  of  the  inverse  problem 
for  this  example.  The  ray  parameter  ranges  from  .12  to  .18  sec/km.  If  a 
sampling  interval  of  &p  =  .01  sec/km  is  used,  the  inverse  problem 


is  set  up  to  estimate  six  values  of  x^  (x(.18)  is  zero): 


X  -  fx>i  9  9  •  •  •  xj 


The  matrix  A  and  vectors  d  and 


— + 

to  are 


(id 


Row 

* 

1 

.01 

.005 

2 

.0? 

.005 

3 

.01 

.01 

.01 

.005 

4 

.01 

.01 

.01 

.005 

5 

.01 

.01 

.01 

.01 

.01 

6 

.01 

.01 

.01 

.01 

.01 

7 

1. 

8 

1. 

9  ' 

1. 

10 

1. 

u 

1. 

12 

13 

-.011 

.075 

14 

s  -  .021 

-.011 

.077 

15 

-  .007 

-.021 

-.012 

.080 

16 

-.004 

-.007 

-.022 

-.01? 

.083 

17 

-.003 

-.004 

-.008 

-.023 

-.01? 

18 

-.252 

.203 

19 

-.020 

-.231 

.185 

20 

.035 

-.018 

-.210 

.168 

21 

.008 

.032 

-.016 

-.189 

.151 

22 

.003 

.007 

.028 

-.014 

-.170 

23 

.029 

.034 

.043 

.068 

.083 

.005 

.005 


.086 


.134 


.15 

.25 

.75 

.85 

1.75 

1.85 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 


-1. 

-T. 

-1. 

-1. 

-1. 

-1. 

-I 

-L 

-L 

-L 

-L 

-L 

-l 

-L 

-L 

-L 

-L 

-.001 
-.001 
-.001 
-.001 
-.001 
-.01  . 


.001 

.001 

.001 

.001 

.00) 

.01 


(12) 
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The  precedence  of  each  equation  is  determined  by  its  corresponding  weighting 
function  u>:.  Here  L  and  e  are  very  large  and  small  numbers,  respec¬ 
tively,  which  serve  to  make  the  corresponding  equations  into  inequalities. 

Rows  1  through  6  impose  the  travel  time  observations.  Notice  that  each  observa¬ 
tion  is  entered  twice  in  order  to  incorporate  the  range  of  uncertainty  in 
the  observation.  Rows  7  through  12  impose  the  condition  _>  0.  Rows  13 
through  17  impose  the  condition  that  z  be  single  valued  in  V,  i.e.,  zi  _>  z^  ^ 
The  coefficients  for  these  rows  were  obtained  by  taking  the  first  differences 
of  the  Weichert-Herglotz  coefficients  defined  in  Equation  (10).  Rows  18  through 
22  suggest  that  in  the  absence  of  other  constraints,  the  velocity  depth  function 
should  be  linear.  These  rows  were  obtained  by  taking  second  differences  of  the 
Weichert-Herglotz  coefficients.  Finally,  the  last  row  requests  that  the 
solution  minimize  the  depth  corresponding  to  p  =  .13  sec/km. 

The  Weichert-Herglotz  weighting  coefficients  for  p  -  .13  sec/km  are 
shown  in  Fig.  1(b).  (A  sampling  interval  of  Ap  =  .005  sec/km  was  used  for 
these  calculations.)  The  minimum  depth  condition  for  p  =  .13  requires 
that  the  weighted  sum  of  the  solution  Xj  over  the  p  interval  [.13,  .18] 
be  zero,  since  the  corresponding  depth  value  in  the  d  vector  has  been 
set  to  zero.  To  obtain  the  velocity  model  which  gives  a  maximum  depth 
for  p  =  .13  sec/km,  the  inverse  problem  is  solved  with  a  very  large 
depth  value  in  the  data  vector  d.  The  p(x)  curves  which  give  the 
maximum  and  minimum  depths  are  shown  in  Fig.  1(b).  In  the  region  around 
p  =  .13  sec/km,  the  solutions  behave  as  would  be  expected,  considering 
the  nature  of  the  extremal  depth  constraints.  At  about  p  =  .15  and 
p  =  .125  sec/km,  the  two  solutions  cross  one  another.  This  is  because  the 
three  travel  times  constrain  the  areas  under  the  p(x)  curves  over  the 
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intervals  [.16,  .18],  [.14,  .18],  [.12,  .18]  and  the  large  excursions 
at  p  =  .13  sec/km  must  be  compensated  elsewhere.  The  velocity-depth 
models  show  a  similar  cross-over  behavior.  The  Solution  which  minimizes 
the  depth  at  p  =  .13  sec/km  has  maximum  depths  at  p  =  .16  and  p  =  .12 
s^c/km.  However,  the  depths  for  ray  parameters  other  than  .13  sec/km  are 
not  necessarily  extremal.  We  have  determined  only  two  of  all  the  possible 
earth  structures  consistent  with  these  data  constraints.  The  v(z) 
envelope  which  encompasses  all  these  models  is  obtained  by  repeating  the 
inverse  problem  with  extremal  depth  conditions  applied  to  each  discrete  p. 

J 

in  turn. 

Note  that  both  v(z)  solutions  in  Fig.  1(c)  have  regions  where  z 
is  constant  in  v.  This  shows  that  the  single-valued  velocity  criterion, 
which  was  imposed  as  an  inequality  equation,  was  actually  solved  as  an 
equality  in  the  final  solutions.  Without  this  constraint,  the  v(z) 
functions  would  not  be  physically  meaningful.  The  single-valued  velocity 
condition  controls  the  shape  of  the  reverse  branches  >  o)  of  the  p(x) 
functions.  The  small  irregularities  in  the  reverse  branches  are  due  to  the 
discrete  nature  of  the  Weichert-Herglotz  integral.  Note  that  some  of  the 
linear  v(z)  relations  are  also  satisfied  in  the  final  solutions.  The 
linear  portions  of  the  v(z)  curves  correspond  to  the  forward  branches 

(Hx  <  the  PM  curves- 

The  travel  time  curves  and  synthetic  seismograms  for  the  two  extremal 
earth  models  are  presented  in  Fig.  1(d).  The  calculated  travel  times  are 
not  consistent  with  the  "observations"  because  we  have  not  used  an  envelope 
to  limit  the  individual  distances  x..  So  far,  we  have  only  placed  integral 
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constraints  on  the  p(x)  curves,  which  is,  in  effect,  the  same  as  using 
three  x(p)  data.  At  those  points  on  the  travel  time  curves  in  Fig.  1(d) 
where  the  slopes  are  .16,  .14  or  .12  sec/km,  the  corresponding  delay  time 
values  are  correct.  However,. the  T(x)  values  are  incorrect  because  these 
slopes  do  not  occur  at  distances  of  20,  40  or  60  km,  as  they  should. 

Fig.  2  shows  the  solutions  which  maximize  the  minimize  the  depths  for 

each  discrete  ray  parameter  between  p  =  .12  and  p  =  .16  sec/km,  with 
Ap  =  .005  sec/km.  The  only  difference  between  each  of  these  solutions  is  in 
the  row  containing  the  Weichert-Herglotz  coefficients  in  the  inverse  problem. 
The  width  of  the  v(z)  envelopes  in  Fig.  2(b)  changes  substantially  with  p. 
This  behavior  is  controlled  by  the  p  limits  on  the  travel  time  integrals. 
The  v (z )  envelope  is  narrower  for  p  values  close  to  these  limits.  Many 

of  the  velocity-depth  models  cross  from  one  extremal  limit  to  the  other  for 

different  values  of  p.  Note,  however,  that  one  particular  solution  cannot 
exhibit  maximum  or  minimum  depths  for  every  p,  since  this  would  violate 
the  travel  time  integrals. 

7  EFFECTS  OF  OBSERVATIONAL  ERRORS  AND  NUMBER  OF  DATA 

Intuitively,  we  expect  the  width  of  the  velocity-depth  envelope  to 
reflect  the  amount  of  uncertainty  associated  with  each  data  point.  Obser¬ 
vation  error  can  be  incorporated  into  the  inverse  problem  in  the  following 
manner:  First,  an  estimate  of  the  uncertainty  +  Ad^  associated  with 
each  data  value  d^  is  made.  For  example,  in  the  previous  problem,  each 
travel  time  observation  was  assigned  an  uncertainty  of  +  .05  sec.  Then, 
two  equations  for  each  data  constraint  are  simultaneously  added  to  the 
inverse  problem;  one  containing  the  value  (d^  +  Ad^)  and  the  other  with 
(d^  -  Ad^)  in  the  data  vector  <J.  The  extremal  depth  constraint  causes 


the  final  error  in  the  L-|  algorithm  to  be  minimized  at  either  one  of 
the  endpoints.  The  algorithm  automatically  selects  the  best  equation 
to  be  exactly  satisfied.  In  general,  the  minimum  depth  solution  selects 
the  one  endpoint  and  the  maximum  depth  solution  the  other.  Hence, 
the  larger  the  uncertainty  in  the  data,  the  greater  the  difference  between 
the  two  solutions.  The  effect  of  data  uncertainty  in  the  previous  problem 
can  be  seen  in  Fig.  2(b),  where  the  v(z)  envelope  for  travel  time  uncer¬ 
tainties  of  +  .1  sec.  is  also  shown.  However,  since  the  widths  of  the 
two  v(z)  envelopes  are  nearly  identical,  we  can  conclude  that  the  large 
depth  limits  in  this  example  are  caused  by  factors  other  than  observation 
error. 

The  next  example  shows  how  the  number  of  observations  affects  the 
resolution  of  earth  structure.  It  is  essentially  a  repeat  of  the  previous 
problem,  except  that  six  travel  time  constraints  are  used.  Their  associated 
p  values  are  .12,  .13,  .14,  .15,  .16  and  .17  sec/km.  Again,  extremal  depthts 
for  each  discrete  ray  parameter  p.  in  the  interval  [.12,  .16]  are  found 

J 

in  turn.  The  solutions  are  presented  in  Fig.  3.  The  range  in  velocity- 
depth  models  is  much  smaller  than  in  the  previous  example.  In  particular, 
the  number  of  "steps"  in  the  v(z)  envelope  has  increased.  The  T(x) 
and  p(x)  solutions  are  also  better  behaved  than  the  ones  in  Fig.  2. 

However,  there  is  still  a  considerable  discrepancy  between  the  observed 
and  calculated  travel  times.  The  next  step  in  the  interpretation  scheme 
is  to  limit  the  p(x)  solutions  in  such  a  way  that  the  corresponding 
velocity-depth  models  become  consistent  with  the  observed  seismograms. 
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8  MATCHING  THE  OBSERVED  SEISMOGRAMS 

In  order  to  depict  accurately  the  range  of  earth  structures  consis¬ 
tent  with  a  given  data  set,  the  inverse  problem  requires  more  than  simple 
linear  constraints  on  the  p{x)  curve.  As  we  saw  in  the  previous  examples, 
the  arrivals  corresponding  to  the  initial  p(x)  solutions  have  travel  times 
and  amplitudes  which  are  radically  different  from  the  observed  seismograms. 
Travel  time  inconsistencies  can  be  eliminated  by  constructing  bounds  in 
the  p-x  plane  and  thereby  restricting  the  range  of  allowable  travel  time 
integrals  for  the  p(x)  solutions.  For  example,  the  travel  time  for  a 
synthetic  arrival  which  occurs  earlier  than  indicated  by  the  observed  data 
can  be  increased  in  the  next  iteration  by  imposing  bounds  which  limit  the 
minimum  values  attained  by  the  p(x)  solutions. 

Amplitude  information  is  also  incorporated  into  the  inverse  problem 
through  the  p(x)  envelope.  The  objective  is  to  make  the  amplitudes  of 
the  synthetic  seismograms  match  the  amplitude  vs  distance  behaviour  of  the 
recorded  data.  The  p(x)  envelope  is  used  to  restrict  the  slopes  of  the 
p(x)  solutions  and  to  control  the  locations  and  lateral  extensions  of  any 
triplications.  Since  recorded  seismograms  usually  exhibit  a  "scatter"  in 
amplitude  between  adjacent  traces,  an  attempt  is  made  to  match  only  those 
amplitude-distance  variations  which  can  be  attributed  to  gross  changes  in 
the  velocity-depth  function. 


18 


«V--  jv.s.* 


The  strategy  for  matching  the  travel  times  and  amplitudes  of  the 
observed  seismograms  is  interactive.  The  first  step  is  to  obtain  a  set 
of  extremal  solutions  with  no  p(x)  envelope  constraint.  We  then  examine 
the  synthetic  seismograms  corresponding  to  each  of  the  solution  models. 

Some  of  the  synthetic  arrivals  will  be  consistent  with  the  observations 
and  some  will  not.  If  we  were  to  color  the  sections  of  the  p(x)  curves 
that  correspond  to  inconsistent  arrivals,  we  would  conclude  that  certain 
portions  of  the  p-x  plane  should  be  excluded  in  order  to  obtain  the 
desired  time-amplitude  behaviour.  If  we  then  incorporate  such  limits 
on  the  p(x)  solutions,  we  can  perform  another  iteration  and  continue 
this  process  until  all  of  the  models  have  a  satisfactory  distribution  of 
arrivals. 

Obviously,  this  approach  is  fairly  crude.  There  may  be  types  of 
amplitude  variations  that  are  not  simply  correlated  with  areas  in  the  p-x 
plane.  The  approach  used  here  could  not  address  such  possibilities.  Our 
experience  from  initial  applications  of  this  technique  is  that  the  use  of 
p(x)  limits  has  handled  the  incorporation  of  amplitudes  to  the  degree  that 
we  are  willing  to  trust  the  amplitude  variations. 

To  illustrate  the  effects  of  the  p(x)  envelope  constraint,  the 
T(x)  data  in  Table  1,  which  were  used  as  data  constraints  for  the  solution 
in  Fig.  2,  have  been  inverted  with  limits  on  the  p(x)  curves.  The  results 
are  presented  in  Fig.  4.  Two  things  are  immediately  obvious:  The  v(z) 
envelope  is  much  narrower  than  the  one  in  Fig.  2(b)  and  the  synthetic 
seismograms  are  now  very  similar  to  the  "observed"  seismograms  in  Fig.  l(z). 
The  spurious  arrivals  have  been  completely  eliminated.  Again,  maximum  and 
minimum  depths  for  p  only  in  the  range  p  =  .12  to  p  =  .16  sec/km  have 
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been  determined.  Also  shown  in  Fig.  4(b)  are  the  extremal  depth  limits  for 
a  slightly  wider  p ( x )  envelope.  There  is  a  close  relationship  between  the 
width  of  the  p(x)  envelope  and  the  range  of  allowable  earth  models. 

9  TRIPLICATION  EXAMPLE 

The  next  example  is  more  realistic.  Fig.  5(a)  shows  four  seismo¬ 
grams  computed  from  the  p(x)  curve  in  Fig.  5(b).  The  earth  model  has  a 
high  velocity  gradient  at  a  depth  of  about  10  km  which  causes  a  triplication 
in  the  p(x)  and  T(x)  curves.  This  example  is  different  from  the  previous 
one  because  all  of  the  arrivals  have  associated  p  values  which  fall  on  the  for¬ 
ward  branches  (dp/dx  <  0)  of  the  p(x)  curve.  These  regions  of  the  p(x)  curve 
are  therefore  quite  well  constrained  while  the  solutions  of  the  range  p  =  .14 
to  .16  sec  km  are  under-constrained.  Fig.  6  shows  the  extremal  depth  solutions 
for  each  discrete  ray  parameter  between  p  =  .12  and  .18  sec/km.  The  uncer¬ 
tainty  in  the  solutions  is  strongly  dependent  on  p.  The  v(z)  envelope 
in  Fig.  6(b)  shows  no  indication  of  the  large  velocity  gradient  at  10  km 
depth.  However,  the  T(x)  curves  in  Fig.  6(c)  do  indicate  the  presence  of 
a  triplication.  Also  shown  in  Fig.  6(b)  is  the  v(z)  envelope  obtained 
when  the  differential  travel  time  constrains  are  omitted.  These  data  have 
a  significant  effect  on  the  width  of  the  velocity-depth  limits. 

Fig.  7  shows  the  extremal  depth  solutions  obtained  using  the  same 
data  as  above  with  an  envelope  limiting  the  p(x)  curves.  The  range  in 
v ( z )  models  is  decreased  considerably  relative  to  the  results  in  Fig.  6(b). 

The  travel  times  and  amplitudes  are  now  quite  similar  to  the  "observations" 
in  Fig.  5(z).  Even  though  the  p(x)  envelope  here  is  quite  narrow,  there 
are  still  some  obvious  discrepancies  between  the  observed  and  calculated 
travel  times  of  the  second  arrivals.  Thus,  the  p(x)  limits  could  be  made 
even  narrower.  Also  presented  in  Fig.  7(b)  is  the  v(z)  envelope  obtained 

by  using  a  slightly  wider  p(x)  envelope. 
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10  OFFSHORE  EXAMPLE 

As  a  final  example,  we  present  the  results  of  our  interpretation 
scheme  applied  to  oceanic  crustal  data  from  Spudich,  et  al.  (1979).  These 
data  have  been  interpreted  by  Spudich,  et  al.  (1979)  by  matching  the 
waveforms  of  the  observed  seismograms  and  synthetics  computed  using  the 
reflectivity  method.  They  have  also  interpreted  these  data  using  the 
travel  time  inversion  scheme  of  Garmany,  et  al.  (1979).  This  technique 
uses  linear  programming  to  invert  for  a  velocity-depth  distribution  directly 
from  a  T(p)  envelope.  Fig.  8  shows  the  final  solutions  obtained  by 
applying  our  interpretation  scheme  to  the  same  x(p)  data  used  by  Spudich, 
et  al.  (1979).  The  p(x)  solutions  have  been  forced  to  lie  within  a  narrow 
envelope  in  order  to  match  the  first  arrival  times  and  amplitudes  of  the 
observed  data.  The  travel  times  for  the  later  arrivals  could  not  be  deter¬ 
mined  from  the  observed  seismograms  with  any  reasonable  accuracy  and  so  are 
not  shown.  Because  these  data  are  lacking,  the  amplitude  information  is 
essential  for  determining  the  nature  of  the  triplications  in  this  example. 

The  relative  amplitudes  shown  in  Fig.  8(d)  are  the  maximum  displacements 
for  each  of  the  observed  seismograms  and  for  synthetics  which  were  computed 
from  the  p(x)  solutions  using  Disc  Ray  Theory  integration. 

Fig.  8(e)  shows  the  velocity  models  for  the  offshore  data  obtained 
by  Spudich,  et  al.  (1979)  and  by  our  interpretation  method.  The  depth  limits 
obtained  using  the  travel  time  inversion  scheme  of  Garmany,  et  al.  (1979)  are 
somewhat  larger  than  those  obtained  using  our  procedure  in  the  first  iteration 
(e.g.,  without  a  p(x)  envelope).  We  have  attributed  this  discrepancy  to 
two  factors;  the  travel  time  inversion  scheme  of  Garmany,  et  al.  (1979)  does 
not  utilize  an  explicit  smoothing  constraint  such  as  a  linear  velocity-depth 
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relationship,  and  the  earth  model  used  by  Spudich,  et  al .  (1979)  for  their 
inverse  problem  had  more  degrees  of  freedom  than  ours  did. 

11  DISCUSSION 

We  have  presented  a  technique  for  interpreting  body-wave  seismograms 
which  involves  modifying  the  p(x)  curve  until  a  best  fit  is  obtained  for 
travel  times  and  amplitudes.  The  p(x)  function  is  represented  by  the 
discrete  distances  x.,  j  =  1,  M  corresponding  to  finite  samplings  of  the 

J 

ray  parameter,  p.,  j  =  1,  M.  Up  to  this  point  we  have  not  mentioned  any 

J 

criterion  for  choosing  a  sampling  interval  Ap.  Ideally,  the  discretization 
of  the  p(x)  curve  should  not  affect  the  solutions.  In  other  words,  the 
sampling  interval  should  be  small  enough  to  make  the  p(x)  curve  appear 
continuous.  In  practice,  the  number  of  degrees  of  freedom  in  the  inverse 
problem  is  limited  by  computer  time.  The  computing  effort  for  the  inverse 

3 

problem  is  proportional  to  NM  ,  where  M  is  the  number  of  unknowns  and  N 
is  the  number  of  rows  in  the  transfer  matrix  A.  An  upper  bound  for  Ap 
for  each  problem  is  determined  by  the  complexity  of  the  p(x)  curve.  The 
linear  interpolation  of  p(x)  is  an  implicit  smoothing  constraint  which 
increases  in  significance  as  the  sampling  interval  is  increased.  In 
general,  Ap  should  vary  as  a  function  of  p.  In  the  examples  presented 
above,  the  sampling  interval  was  deliberately  kept  large  in  order  to  keep 
the  diagrams  readable. 

We  feel  that  there  are  two  aspects  of  this  interpretation  procedure 
which  require  additional  research.  The  first  involves  extending  the  scheme 
to  include  earth  structures  containing  low-velocity  zones.  The  second  is 
related  to  the  method  for  determining  the  p(x)  envelope,  which,  at  the 
present  time,  requires  human  intervention.  We  found  that  one  quickly 
learns  how  to  manipulate  the  p(x)  curve  in  order  to  effect  desired  changes 


22 


in  the  travel  times  and  amplitudes  of  the  synthetic  seismograms.  However, 
with  some  additional  programing,  this  step  could  be  made  automatic. 
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TABLE  CAPTIONS 

Table  1.  Travel  time  data  for  the  three  seismograms  in  Fig.  1(a). 


FIGURE  CAPTIONS 


Figure  1.  Two  solutions  to  the  inverse  problem  which  give  extremal 

depths  for  one  velocity,  (a)  shows  three  seismograms  which 
provide  the  data  for  the  inversion.  The  solid  circles  represent 
travel  time  uncertainties  of  +  .05  sec.  The  seismograms  were 
computed  from  the  monotonic  p(x)  function  in  (b)  (solid  curve 
which  has  a  corresponding  earth  model  depicted  as  a  solid  curve 
in  (c).  The  Weichert-Herglotz  coefficients  for  p  =  .13  used 
in  the  extremal  depth  condition  are  shown  schematically  as  solid 
circles  in  (b).  The  short  and  long  dashed  curves  in  (b)  are 
the  p(x)  solutions  which  minimize  and  maximize  respectively 
the  depth  for  p  =  .13  sec/km.  (c)  shows  the  corresponding 
velocity-depth  models.  The  T(x)  curves  and  synthetic  seismo¬ 
grams  for  the  extremal  depth  solutions  are  presented  in  (d). 

Note  that  that  travel  times  and  amplitudes  for  these  solutions 
are  not  consistent  with  the  data. 

Figure  2.  Extremal  depth  solutions  for  each  discrete  ray 

parameter  between  .12  and  .16  sec/km.  The  data  consist  of 
three  travel  times  with  uncertainties  of  +  .05  sec,  as  for 
the  results  in  Fig.  1.  The  dashed  curves  in  (b)  represent 
the  depth  limits  obtained  by  using  three  travel  time  data 
with  uncertainties  of  +  .1  sec.  Note  the  discrepancies 
between  the  T(x)  curves  in  (c)  and  the  data  in  Fig.  1(a). 


Figure  3.  Extremal  depth  solutions  for  six  travel  time  data  with 


uncertainties  of  +  .05  sec.  The  "true"  earth  model  for 
this  example  is  the  solid  v(z)  curve  in  Fig.  1(c).  The 
only  difference  between  the  solutions  presented  here  and  those 
in  Fig.  2  is  in  the  number  of  travel  time  data  used  to  constrain 
the  p(x)  functions.  Note  that  the  extremal  depths  in  (b)  are 
much  narrower  than  the  ones  obtained  using  three  T(x)  data, 
presented  in  Fig.  2(b).  The  dashed  curves  in  (b)  represent  the 
depth  limits  obtained  by  using  travel  time  uncertainties  of 
+  .1  sec. 

Figure  4.  Extremal  depth  solutions  for  the  three  travel  time  data  in 

Fig.  1(a)  with  a  p(x)  envelope  constraint  added.  The  T(x) 
data  have  uncertainties  of  +  .05  sec.  The  v(z)  envelope  in 
(b)  is  much  narrower  than  the  one  obtained  without  using  a  p(x) 
envelope,  presented  in  Fig.  2(b)  .  The  travel  times  and  synthetic 
seismograms  in  (c)  are  quite  similar  to  the  data  in  Fig.  1(a). 

The  dashed  curves  in  (b)  are  the  depth  limits  obtained  by  using 
the  p(x)  envelope  shown  as  a  dashed  curve  in  (a). 

Figure  5.  (a)  shows  four  seismograms  which  provide  the  data  for  the 

"triplication"  example.  The  "true"  p(x)  and  v(z)  functions 
are  presented  in  (b)  and  (c).  The  travel  time  uncertainties 
are  +  .05  sec. 


Figure  6.  Extremal  depth  solutions  for  the  travel  time  data  in  Fig.  5(a). 

No  envelope  was  used  to  restrict  the  p(x)  curves.  The  data 


constraints  for  these  results  consist  of  four  first  arrival 
times  and  two  differential  travel  times.  The  depth  limits 
shown  as  dashed  curves  in  (b)  were  obtained  using  the  four 
first  arrival  time  constraints  only. 

Figure  7.  Extremal  depth  solutions  for  the  data  in  Fig.  5(a)  with  a  p(x) 
envelope  constraint  added.  Note  how  narrow  the  depth  limits 
in  (b)  are  compared  to  those  obtained  without  using  a  p(x) 
envelope  (Fig.  6(b)).  The  first  arrival  times  for  the  above 
solutions  are  comparable  to  the  "observed"  values  in  Fig.  5(a) 
However,  the  p(x)  envelope  must  be  made  even  narrower  in 
order  to  fit  the  travel  times  for  the  later  arrivals.  The 
dashed  curves  in  (b)  represent  the  depth  limits  obtained  by 
using  the  dashed  p(x)  envelope  in  (a). 

Figure  8.  Results  obtained  by  inverting  oceanic  c.ustal  data.  The 

v(z),  p(x),  T(x)  and  amplitude  functions  obtained  in  the  final 
iteration  are  presented  as  solid  curves  in  (a),  (b),  (c)  and 
(d)  respectively.  The  final  p(x)  envelope  for  these  solutions 
can  be  seen  in  (b).  The  first  depth  to  be  maximized  occurs  at 
v  =  4  km/sec  and  the  next  one  is  at  5.9  km/sec.  Between  5.9 
and  8  km/sec  the  velocity  increments  are  quite  small.  The 
solid  circles  in  (c)  represent  the  errors  for  the  observed 


(Figure  8  continued) 


travel  times.  The  maximum  amplitudes  for  the  observed 
seismograms  are  depicted  as  solid  circles  in  (d).  (e)  shows 

the  depth  limits  obtained  by  the  travel  time  inversion  scheme 
of  Garmany,  et  al.  (1979)  ("+"'s)  and  by  our  technique,  with 

and  without  using  the  p(x)  envelope  constraint  (long  dashed  and 
solid  curves  respectively).  The  layered  crustal  model  obtained 
by  Spudich,  et  al.  (1979)  by  waveform  fitting  is  shown  in  (e) 
as  a  short  dashed  curve. 
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